##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)
library("scales")

##########################################################################################

option_list <- list(
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--molecular_type"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){

    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    info_file <- paste(work_dir,"/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv",sep="")
    images_path <- paste(work_dir,"/images/mutBurden",sep="")
    molecular_type <- "All"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

molecular_type <- opt$molecular_type
info_file <- opt$info_file
images_path <- opt$images_path

dir.create( images_path , recursive = T )

###########################################################################################

dat_sample <- data.frame(fread( info_file ))

###########################################################################################
if(molecular_type=="All"){
    mss_sample <- subset(dat_sample , TCGA_Class %in% c("GS" , "CIN"))$Patient
    dat_sample <- subset( dat_sample , Patient %in% mss_sample )
}else if(molecular_type=="GS"){
    gs_sample <- subset(dat_sample , TCGA_Class=="GS")$Patient
    dat_sample <- subset( dat_sample , Patient %in% gs_sample )
}else if(molecular_type=="CIN"){
    cin_sample <- subset(dat_sample , TCGA_Class=="CIN")$Patient
    dat_sample <- subset( dat_sample , Patient %in% cin_sample )
}

## 多个样本负荷取中位数
dat_plot2 <- dat_sample %>%
group_by( Patient , Class , Type , Tobacco , Alcohol , HP ) %>%
summarize( BurdenAll = median(BurdenAll) ,BurdenExon = median(BurdenExon) )

###########################################################################################

base_col <- c("Tobacco" , "Alcohol" , "HP")
dat_plot <- subset(dat_plot2 , Type != "IM + IGC + DGC" )
dat_plot$Tobacco <- factor( dat_plot$Tobacco , levels = c("No" , "Smoke") )
dat_plot$Alcohol <- factor( dat_plot$Alcohol , levels = c( "No" , "Drink") )
dat_plot$HP <- factor( dat_plot$HP , levels = c("Negative" , "Positive") )

result <- c()
for( base_use in base_col ){
    ## 单因素
    # logistic回归模型拟合
    print(base_use)
    dat_plot$use_col <- dat_plot[[base_use]]
    dat_plot$useCol <- dat_plot[["BurdenExon"]]
    
    ## 分不同时期计算
    for( class in unique(dat_plot$Class) ){
        tmp_dat <- subset( dat_plot , Class == class )
        tmp_dat <- tmp_dat[!is.na(tmp_dat$use_col),]
        tmp_dat <- subset( tmp_dat , use_col!="NA" )
        model <- glm( useCol ~ use_col , data = tmp_dat)

        ## p值
        P <- summary(model)$coef[2,4]
        ## 95%可惜区间
        OR <- cbind(OR = coef(model), confint(model))[2,]
        tmp <- data.frame( baseType = base_use , model = "single" , var = gsub( "use_col" , "" , rownames(summary(model)$coef)[2]) , time = class , p = P , coef = OR[1] , coef_025 = OR[2] , coef_975 = OR[3] )

        result <- rbind(result , tmp)
    }
}


out_name <- paste0( images_path , "/baseInfo_mutBurden.",molecular_type,".tsv" )
write.table( result , out_name , row.names = F , quote = F , sep = "\t" )
